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1)  Objectives 


The  primary  objective  of  the  computational  project  is  to  develop  a  computer 
model  to  simulate  the  continuous  production  of  a  particulate  solid  film  from  a  low 
density  feed  phase,  by  uniaxial  compaction. 

In  the  common  process  for  the  continuous  production  of  a  binary  composite 
material,  a  powder  is  fed  from  a  hopper,  a  second  powder  is  fed  from  another  hopper, 
both  may  mixed  and  with  a  liquid  feed,  or  dry  with  gravity  or  a  gas  cyclone.  The  mixed 
colloidal  suspension  is  then  piped  for  processing  to  a  solid  material.  Typically,  at  some 
later  stage,  the  particles  in  the  binary  suspension  are  to  be  compacted  to  a  dense  solid 
state.  The  imiaxial  compaction  process  could  be  a  sedimentation,  centrifugation, 
filtration,  or  slip  casting  as  in  ceramics,  or  tabletting  as  in  pharmaceuticals,  or  pressure 
injection  molding.  Other  processes  for  producing  thin  solid  films  involve  deposition  in 
large  thermal  gradients  and  splat  quenching  on  to  a  colder  substrate. 


Figure  1  :  Schematic  chart  of  a  typical  production  process  for  a  binary  composite  material 
processed  continuously  via  a  colloidal  suspension. 
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The  shear  flow  and  deformation  processes  that  determine  the  final  product  require  a 
prediction  of  the  fluid  mechanics  of  the  mixed  suspension  [2]  The  rate  of  flow  and  the 
distribution  of  particulates  as  a  function  of  flow  rate,  i.e.  the  fluid  mechanics  of  the  colloidal 
suspension  in  a  given  geometry  is  crucially  dependent  of  the  solid  fractions  of  the  two  powders 
and  their  variable  particle  size  distributions,  and  on  the  viscosity  and  temperature  of  the  liquid. 
Control  of  the  composition  of  the  feeds  and  the  liquid  flow  therefore  requires  feedback 
information  of  the  constitutive  rheological  relations  that  will  enable  computational  fluid 
mechanics  of  the  suspension  after  mixing  and  during  processing. 

When  a  colloidal  suspension  is  subjected  to  any  stress  deformation,  the  osmotic 
thermodynamic  and  transport  properties  of  the  particles  themselves  become  a  function  of 
deformation  rate,  in  addition  to  the  stress.  Finite  element  and  finite  difference  CFD  simulations 
of  particulate  systems  require  the  osmotic  pressure,  the  particle  diffusivity,  and  the  particle 
kinetic  energy  conductivity,  all  as  a  function  of  the  rate  of  deformation. 

In  principle  there  are  three  possible  ways  to  obtain  all  or  part  of  this  information.  The 
conventional  chemical  engineering  approach,  either  batch-sample  or  on-line,  is  to  measure  the 
“flow  curve”  of  the  suspension,  i.e.  the  stress  as  a  function  of  concentrations,  and  the  rate  of 
shear  strain,  concentrations  and  temperature,  and  to  utilize  that  information.  This  approach  is 
intrusive  and  fraught  with  problems.  Under  shear  flow  the  particles  reorganize  themselves  to 
minimize  the  stresses.  Rheometric  viscosity  data  for  dense  suspensions  are  geometry  dependent, 
and  are  neither  scalable  nor  transferable  from  one  rheometric  or  processing  device  to  another 
[2,3]  Elongational  viscosities  are  impossible  to  access  experimentally. 

An  alternative  approach  is  to  obtain  all  the  requisite  details  of  the  particle  size  and  shape 
distributions,  the  liquid  medium  density  and  viscosity,  etc.,  and  to  utilize  the  methods  of 
computational  statistical  mechanics.  In  principle,  given  the  essential  information,  mesoscale 
modeling  will  eventually  be  able  to  make  ab  initio  predictions  of  the  fluid  mechanics  of  complex 
suspensions,  but  there  is  a  long  way  to  go.  The  programs,  which  have  been  developed  here  for 
compaction  processes,  can  be  regarded  as  a  preliminary  step  along  the  path  to  an  integrated 
computational  rheometry  of  colloidal  systems. 

In  recent  years  there  have  been  developments  in  non-intrusive  methods  that  gather  data 
by  various  means,  and  utilize  that  advances  in  computing  to  reconstruct  a  detailed  picture  of  the 
interior.  These  methods  are  also  promising  to  obtain  information  on  the  concentration  of 
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particles  and  velocities  of  particles  during  shear  flow  and  processing  and,  potentially,  can  yield 
the  necessary  constitutive  rheological  relations,  by  on  line  techniques.  One  such  method  is  the 
use  of  Raman  spectra  tomography,  which  is  being  pioneered  in  the  Materials  Directorate  AFRL 
in  Dayton.  A  second  and  subsidiary  objective  in  this  project  is  to  test  a  literature  method  for  the 
production  of  a  prototype  monodisperse  latex,  with  the  particle  diameter  of  the  wavelength  of 
green  light  (532nm),  in  laboratory  experiments. 

2)  Mesoscale  model 

The  simulation  of  colloidal  particle  motions  by  methods  akin  to  Monte  Carlo  and 
molecular  dynamics  of  computational  statistical  mechanics  have  come  to  be  known  as  mesoscale 
modeling  [4].  Both  granular  dynamics  simulations  of  inelastic  particles,  and  Stokesian  dynamics 
simulations  of  particles  in  hydrodynamic  media  are  more  akin  to  NEMD  (non-equilibrium 
molecular  dynamics)  in  that  there  is  no  inherent  thermodynamic  equilibrium.  In  the  simulation 
of  flows  and  deformations  of  colloidal  systems,  the  particle  movements  are  driven  by  externally 
imposed  forces. 

If  we  have  sufficient  knowledge  of  the  interparticle  forces,  and  other  essential  particle 
characteristics,  such  as  polydispersity  of  size  and  shape,  surface  charge  and  friction  effects  etc. 
in  principle  the  entire  constitutive  rheological  relations  that  are  required  for  fluid  mechanics 
predictions  in  given  geometries  can  be  obtained.  “Computational  rheometry”  is  still  in  very  early 
stages.  In  principle  it  is  possible  because  the  two  major  obstacles  to  obtaining  experimental  flow 
curves  via  conventional  rheometry  can  be  eliminated  in  idealized  computer  experiments. 

The  general  methodology  that  has  been  employed  is  classical  molecular  dynamics  in 
which  the  equations-of-motion  of  large  numbers  of  particles  are  solved  simultaneously.  These 
could  be  macromolecules,  nano-clusters,  colloidal  particles  in  suspension,  or  inelastic  granular 
particles  which  dissipate  energy  on  collision.  The  interparticle  forces  vary  with  model  material. 
The  equations  of  motion  can  be  readily  changed  from  Newtonian,  to  granular  with  inelastic 
interactions,  or  Stokesian/Brownian  for  suspensions. 

In  the  present  preliminary  investigation  of  continuous  thin-film  solidification,  the 
equations  of  motion  have  been  so  far  restricted  to  classical  |Newtonian. 
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Sliding-brick  periodic  boundaries  [5] 


The  computational  tricks  are  first  to  eliminate  the  presence  of  driving  surfaces  by 
periodic  boundary  conditions  and  then,  secondly,  to  impose  a  uniform  homogeneous  shear  field 
by  velocity  scaling  methods.  Both  these  conditions  enable  the  simulation  of  truly  homogeneous 
shear  flow  regimes  with  the  system  steady-state  variables  prescribed.  It  should  be  noted  the 
homogeneous  shear  steady-state  computations  have  no  real  experimental  counterpart.  In 
principle,  however,  these  simulations  will  lead  to  the  computation  of  the  constitutive  rheological 
relations  for  the  model  material.  This  information,  in  turn,  will  enable  CFD  (computational  fluid 
dynamics)  predictions  for  colloidal  fluids  in  processing  geometries. 


Figure  2  :  Periodic  boundaries  for  simulation  of  steady-state  homogeneous  flow  of  a  particulate 
suspension;  both  the  linear  velocity  gradient  and  the  concentration  remain  uniform. 

Uniaxial  compaction  periodic  boundaries  [6] 

Periodic  boundary  conditions  have  also  been  devised  for  continuous  elongational  and 
compressive  deformations  and  these  enable  the  computation  of  non-linear  elongational 
properties  of  particulate  systems.  These  boimdary  conditions  can  be  generalized  to  a  whole  range 
of  uniaxial  compaction  processes  and  in  principle  can  be  use  to  simulate  the  processes  referred 
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to  in  the  introduction  for  the  production  of  binary  composite  solids  for  example.  Hence,  the 
constitutive  rheological  relations  for  the  solidification  processes  for  colloids  can,  in  principle,  be 
predicted  computationally. 

boundary  driving 

velocity  (vb)  force  (Fa) 


Figure  3  :  Periodic  uniaxial  compaction  boundaries;  under  certain  conditions  of  the  prescribed 
driving  force  and  interface  velocity,  this  system  exhibits  two-phase  behavior  as  described  in  the 
text.  The  pressures  and  diffusivities  of  the  “feed”  and  the  “bed”  phases,  as  a  fimction  of  the 
driving  force  (Fa)  and  the  rate  constant  (Vb),  are  required  for  CFD  simulations  of  the 
corresponding  batch  processes. 

When  an  external  force,  such  as  a  gravitational  field,  is  applied  to  a  dynamical  system,  if 
the  boundary  in  the  direction  of  the  force  is  rigid,  the  particles  will  condense  and  a  time- 
dependent  process  will  lead  to  a  batch  solid  film  at  that  boundary.  If,  on  the  other  hand,  the 
boundary  is  periodic,  the  particles  will  eventually  just  move  in  a  steady  state  through  the  system, 
fi'om  one  side  to  the  other,  recycling,  with  no  change  in  properties.  When  a  semi-permeable 
hybrid  boimdary  is  used  however,  a  steady-state  uniaxial  compaction  process  is  effected  in  which 
a  low  density  "feed"  phase  moving  with  a  high  mass  velocity  moves  continuously  in  coexistence 
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with  a  high-density  solid  material  phase  moving  with  a  lower  mass  velocity  (to  maintain  the 
mass  and  force  balance),  and  an  interface  between  them. 

Processes  that  do  not  characterize  the  solid  state  generally  produce  non-crystalline  solid 
materials.  For  example,  solids  produced  by  industrial  crystallization  or  precipitation  are 
processed  by  uniaxial  compactions,  such  as  sedimentation,  filtration,  or  centrifugation.  Because 
of  the  change  in  conditions  at  the  deposition  interface  with  time  in  these  batch  processes  (the 
driving  force  usually  weakens  with  time),  an  inhomogeneous  bulk  solid-state  material  is 
produced  that  cannot  be  characterized.  This  implies  that  problems  will  be  encountered  in  defect 
and  grain  boundary  inhomogeneities,  non-reproducibility  and  quality  control. 

These  imperfections  cannot  be  tolerated  in  many  modem  materials.  If  it  were  somehow 
possible  to  devise  well-characterized  steady-state  processes  that  continuously  produce  advanced 
solid  composite  materials  such  as  carbon  composites,  ceramics  etc.,  it  might  eventually  be 
possible  to  uniquely  characterize  these  advanced  solid  materials  by  just  one  (if  gravity  is  fixed) 
or  possible  two  additional  state  variables  that  describe  the  process  that  produces  the  non¬ 
equilibrium  solid. 

The  imique  characterization  of  non-equilibrium  Solid  materials  by  the  thermodynamic 
conditions  at  which  they  exist  (temperature  and  either  pressure  or  density),  and  the  process  that 
produces  them  puts  the  characterization  of  the  state  of  advanced  materials  a  sounder  scientific 
base.  Equation-of-state,  transport  coefficients,  and  other  constitutive  properties  could  be  raised 
to  a  status  of  reproducibility  to  that  of  equilibrium  fluids  and  some  equilibrium  well- 
characterized  crystalline  materials.  A  complete  control  model  of  these  idealized  material 
production  processes  is  then  possible.  In  the  following  section,  the  idealized  process  of  uniaxial 
compaction  is  investigated  in  a  system  of  steady-state  conditions  using  simulation  programs  that 
have  been  developed  as  a  part  of  this  project. 

In  the  case  of  solids  produced  from  granular  materials,  powders,  suspensions  of  powders 
and  other  colloidal  particles,  the  state  of  the  feed  materials  must  also  be  addressed,  as  colloids 
are  generally  difficult  to  characterize.  Below  we  develop  a  program  for  simulating  the 
continuous  production  of  a  high  density  solid  phase  continuously  in  steady-state  equilibrium  with 
a  low  density  feed  phase.  The  determination  of  the  properties  of  these  materials  by 
computational  statistical  mechanics  methods  at  the  colloidal  particulate  level  is  the  immediate 
priority  of  the  present  project. 
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3)  Simulations 


In  all  the  programs  used  and  developed  for  this  project  to  date,  (see  APPENDIX  1)  the  6N  coupled 
first-order  differential  equations  of  classical  particle  motion  are 

dri/dt  =  Vi 

and  dvi/dt  =  Fi  /m 

Simple  one-step  predictor  methods  are  adequate.  In  a  finite  difference  form 

d^rj/dt^  =  Fi(t)  /mi  ={[ri(t+dt)  -  ri(t)]  -  [  ri(t)  -  ri  (t -dt)]}/dt^  +  Odt^ 
where  dt  denotes  the  time  step  in  the  numerical  integration.  The  simple  algorithm  is  then 
ri(t+dt)  =  2  ri(t)  -ri  (t  -  dt)  +  dt^  Fi  (t)  /  mi  +  Odt'' 

The  algorithm  can  be  made  self-starting  if  the  initial  velocities  are  specified  and  used  to  define 
ti  (t  —  dt)  and  specify  at  t=0  i.e.  ri  (t  -  dt)  =  ri(t)  -  Vi  (dt=0)  /dt 
then  the  algorithm  for  position  and  velocities  respectively  are 
ri(t+dt)  =  ri(t)  +  dt  Vi  (t)  +  dt^  Fi  (t)  /  mi  +  Odt** 

Vi(t-i-dt)  =  Vi(t)  +  dtFi(t)/mi 

The  whole  system  is  divided  up  into  cells,  nix,  nly  and  nlz  being  the  total  number  of 
neighborhood  compartments  in  each  dimension.  The  total  number  of  cells  is  NLX  x  NLY  x 
NLZ.  The  particles  are  then  indexed  in  order  of  position  beginning  with  the  first  particle  in  cell 
1,  and  keeping  a  track  of  the  number  of  particles  in  each  cell.  As  the  particles  move  aroimd  the 
system,  their  indices  change.  The  code  for  updating  the  indices  is  much  faster  than  conventional 
LINK  cell  methods  and  may  be  more  amenable  to  parallelization,  (see  APPENDIX  1) 

a)  MNCELL.f 

A  computer  program  has  heen  developed  as  an  extension  of  the  previous  program 
(referred  to  as  LJLINK.f),  for  a  system  of  three  dimensional  particles  that  was  used  to  explore 
the  phase  behavior  and  continuous  production  of  an  amorphous  solid  of  n=12  soft  spheres  by 
Warr  and  Woodcock  [6].  The  program  requires  3  basic  ingredients  over  and  above  a  classical  3D 
homogeneous  MD  program  [4] 
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System  (model)  :  interparticle  forces  for  colloidal  particles  in  the  first  approximation  have 
two  components,  a  steep  positive  repulsion  clouds  overlap  at  short-range,  and  a  longer  range 
negative  attraction  which  decrease  with  an  inverse  power  (1-6)  of  the  distance.  Here  a  more 
general  form  of  the  Lennard-Jones  pair  potential,  m-n  potential,  in  which  the  minimum  distance 
(ro)  is  used,  rather  than  the  collision  distance  of  zero  potential. 

Eq  ~  Enon-bonded  particleas  ^AB  {  (^oAB  /j")  ■2(  roAB  /r)  } 

The  m-n  pair  potential  contains  only  two-parameters  the  energy  at  the  minimum  (-s)  and  the 
minimum  energy  distance  (ro)  when  the  force  is  zero:  n  determines  the  “hardness”  of  the  potential  and 
m  determines  the  “range”  of  the  attraction. 

For  binary  systems  the  Lorentz-Berthelot  rules  for  determining  the  parameters,  that  were 
originally  applied  to  mixtures  of  rare  gas  atoms  are  used 

Sab  ~  (^AA^BB  )  and  Toab  ~  (foAA  Tobb  ) 

For  ionic  systems  Coulombic  potentials  the  charge  distribution  can  be  represented  by  a  set  of 
two  or  more  point  charges,  usually,  but  not  necessarily,  centered  on  the  particle,  but  in  any  case  fixed 
relative  to  the  particle  in  space.  The  total  Coulombic  potential  is  then 

Ec  =  Sij  charge  pairs  ZjZje^j  where  Zi  and  Zj  are  the  charges  in  dimensions  of  the  charge  of  an 
electron  (e).  The  programs  develoed  are  general  but  all  the  test  runs  to  date  refer  to  an  m-n  single 
component  fluid  n=12  and  m=6. 

System-state  variables:  cell  dimensions,  density,  temperature 

For  a  fixlly  periodic  system,  once  the  size  and  shape  is  specified  (  in  the  three  side  length 
variables  XF  YF  ZF  with  XF  =1  the  dimension  of  length  used  in  the  program. 

Volume  =  XF  x  YF  x  ZF 

The  number  of  particles  N  is  fixed 

Mean  density  rho=  NA^ 

Kinetic  energy  =  Zmv^ 

External  driving  force  =  Fd 
Boimdary  velocity  =  Vb 

All  the  variables  are  in  m-n  pair  potential  dimensionless  reduced  xmits. 

It  is  the  two  steady-state  variables  that  conspire  to  create  in  certain  regions  of  parameter 
space  a  twO-phase  solidifying  system  of  a  low  density  phase  moving  with  a  high  velocity  (feed 


10 


phase)  in  coexistence  with  a  high  density  solid  phase  moving  with  a  low  velocity  with  and 
interface. 

Boundary  conditions: 

For  the  test  program  MNCELL  there  is  full  periodicity  in  both  of  the  directions 
perpendicular  to  flow,  i.e.  X  periodic  Y  periodic  Z  semi-permeable  (finite  Vh) 

A  sequence  of  test  runs  on  the  USAFR  Computational  facility  Orion 
(acomputerforme.com)  at  San  Antonio  have  been  earned  out.  The  full  Output  are  on  file  as 
labeled  in  table  1  below. 


TABLE  1 :  Test  runs  for  MNCELL  rho=Nro^A^=  1 .0000/  T  =1 .0000;  cubic 


Run/OP 

N 

Time/dt 

Fd 

Vb 

<U> 

<p> 

Comments 

MNOPl 

2048 

100000 

10.0 

0.1 

-2.277 

30.21 

2-phase/gas-cryst/sharp 

MNOP2 

2048 

40000 

10.0 

0.01 

-2.330 

31.04 

perfect  crystal/zero  flux 

MNOP3 

6912 

100000 

1.0 

1.0 

-4.019 

2.67 

2-phase/  dense  fl.  am.. sol 

MNOP4 

8788 

100000 

10.0 

0.05 

-0.878 

38.21 

Equilibration  run/slow 

MNOP5 

8788 

1000000 

10.0 

0.05 

-0.716 

43.36 

2-phase/gas-cryst/ sharp 

MNOP6 

8788 

1000000 

1.0 

0.05 

-4.488 

6.173 

2-phase/gas-cryst/  diffuse 

MNOP7 

8788 

1000000 

1.0 

0.05 

-4.493 

6.143 

2-phase/gas-cryst/sharper 

MNOP8 

8788 

1 

1000000 

1.0 

0.05 

-4.286 

6.235 

2-phase/gas-cryst/sharp 

At  high  velocities  (>~1)  and  low  driving  force  (<~1)  a  single  amorphous  phase  is 
obtained,  but  the  table  above  shows  that  there  is  a  wide  range  of  parameter  space  that  leads  to  a 
continuous  solidification,  with  a  low  density  feed,  an  interface,  and  a  high  density  solid.  The 
solid  may  be  crystalline,  amorphous  or  possible  metastable. 

Run  3  (Vb=l,Fd=l)  is  2-phase  but  close  to  the  one  phase  boundary. 

For  all  the  runs  profiles  in  the  heterogeneous  Z-direction  (i.e.  the  direction  of  flow)  are 
accumulated.  A  typical  set  of  profiles  for  a  steady-state  system  showing  two  phases  with 
crystallization  is  given  in  Table  2  below. 

TABLE  2:  EXAMPLE  of  2-phase  PROFILES  for  TEST  RUN  1 
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The  profiles  are.  Respectively,  index  of  bin,  height  Z  (if  force  is  gravity),  density(rho), 
velocity  in  z-direction  (of  flow),  mass  flow  rate  (flux),  kinetic  energy  (KE),  potential  energy 


(PE),  and  interfacial  stress  for  each  layer  (-Pxy)z. 


h 

Z/ro 

RHO 

VZ 

<f lux> 

KE 

PE 

STRESS 

50 

12.45 

0.013 

-1.811 

-0.024 

2.379 

-3.758 

-0.087 

49 

12.19 

0.013 

-1.947 

-0.025 

2.702 

-4.416 

-0.024 

48 

11.94 

0.017 

-1.530 

-0.027 

2.066 

-4.013 

-0.011 

47 

11.68 

0.025 

-1.097 

-0.028 

1.607 

-3.439 

-0.050 

46 

11.43 

0.028 

-1.043 

-0.029 

1.658 

-2.800 

-0.067 

45 

11.18 

0.024 

-r.283 

-0.030 

1.927 

-2.621 

-0.033 

44 

10.92 

0.023 

-1.377 

-0.032 

2.000 

-2.675 

-0.018 

43 

10.67 

0.023 

-1.420 

-0.033 

2.019 

-2.736 

-0.022 

42 

10.41 

0.023 

-1.457 

-0.034 

2.072 

-2.773 

-0.016 

41 

10.16 

0.025 

-1.441 

-0.035 

2.139 

-2 . 891 

-0.015 

40 

9.91 

0.025 

-1.460 

-0.037 

2.178 

-3.004 

-0.010 

39 

9.65 

0.026 

-1.472 

-0.038 

2.210 

-3.097 

-0.011 

38 

9.40 

0.027 

-1.469 

-0.039 

2.195 

-3.213 

-0.010 

37 

9.14 

0.027 

-1.503 

-0.041 

2.194 

-3.328 

-0.012 

36 

8.89 

0.028 

-1.480 

-0.042 

2.131 

-3.403 

-0.008 

35 

8.64 

0.029 

-1.475 

-0.043 

2.145 

-3 .476 

-0.014 

34 

8.38 

0.030 

-1.499 

-0 . 044 

2.209 

-3.554 

-0.011 

33 

8.13 

0.031 

-1.484 

-0.045 

2 . 143 

-3.796 

-0.011 

32 

7 . 87 

0.031 

-1.488 

-0.047 

2 . 157 

-3.862 

-0.011 

31 

7.62 

0.034 

-1.430 

-0.048 

2.104 

-3.907 

-0.002 

30 

7.37 

0.038 

-1.308 

-0.049 

2.010 

-3.842 

0.007 

29 

7.11 

0.068 

-0.738 

-0.050 

1.597 

-3.316 

0.020 

28 

6.86 

0.410 

-0.118 

-0.048 

1.106 

-3.137 

0.086 

27 

6.60 

1.539 

-0.033 

-0.051 

1.020 

-3.605 

3.929 

26 

6.35 

0.518 

-0.105 

-0.054 

1.042 

-4.013 

0.692 

25 

6.10 

0.865 

-0.054 

-0.047 

1.021 

-4.986 

-2.228 

24 

5.84 

2 . 896 

-0.017 

-0.050 

0.991 

-5.145 

-1.371 

23 

5.59 

0.528 

-0.106 

-0.056 

0.995 

-5.176 

0.419 

22 

5.33 

0.723 

-0.062 

-0.044 

0.994 

-5.838 

-1.736 

21 

5.08 

3.636 

-0.014 

-0.052 

0.982 

-5.363 

-14 . 090 

20 

4.83 

0.325 

-0.155 

-0.050 

0. 927 

-5.298 

-0.201 

19 

4.57 

0.842 

-0.046 

-0.039 

0.967 

-5.845 

-1.149 

18 

4.32 

3.701 

-0.015 

-0.056 

0.977 

-5.169 

-17.902 

17 

4.06 

0.273 

-0.163 

-0.044 

0.852 

-5.209 

0.160 

16 

3.81 

2.319 

-0.005 

-0.012 

0.977 

-5.073 

-10.525 

15 

3.56 

2.260 

-0.034 

-0.077 

0.982 

-4.674 

-11.028 

14 

3.30 

0.320 

-0.131 

-0.042 

0.837 

-4.523 

-0.063 

13 

3.05 

4.096 

-0.010 

-0.039 

1.025 

-4.690 

-16.131 

12 

2.79 

0.252 

-0.168 

-0.042 

0.843 

-4.258 

-0.375 

11 

2.54 

0 . 823 

-0.025 

-0.020 

0.935 

-4.017 

-1.515 

10 

2.29 

3 . 849 

-0.015 

-0.057 

0.980 

-4.034 

-4.991 

9 

2.03 

0.252 

-0.143 

-0.036 

0.827 

-3.747 

0.447 

8 

1.78 

4.096 

-0.008 

-0.031 

1.042 

-3.509 

3.837 

7 

1.52 

0.230 

-0.169 

-0.039 

0.951 

-3.236 

-0.787 

6 

1.27 

0.879 

-0.010 

-0.009 

0.952 

-2.440 

1.998 

5 

1.02 

3 . 828 

-0.015 

-0.057 

0.987 

-2.354 

12.101 

4 

0.76 

0.134 

-0.225 

-0.030 

1.006 

-4.123 

-0.495 

3 

0.51 

4.096 

-0.007 

-0.030 

1.091 

-1.324 

32.627 

2 

0.25 

0.059 

-0.473 

-0.028 

1.168 

-2.813 

-0.365 

1 

0.00 

4.096 

0.112 

0.458 

1.053 

-1.659 

-25.570 
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b)  MNSLAB.f  (see  appendix  I) 

A  derivative  computer  program  has  been  developed  in  which  the  periodic  boundaries  in 
the  Y-direction  are  made  to  be  rigidly  reflecting  rather  than  periodic.  The  program  the  simulates 
the  continuous  production  of  a  thin  solid  film  in  which  the  particles  are  constrained  to  move 
between  the  two  boundaries. 

The  system  and  system-state  parameters  are  the  same  as  MNCELL  above,  but  the  system 
is  not  periodic  in  one  of  the  perpendicular  directions  (Y).  The  periodic  boundary  conditions  in 
the  Y-direction  are  replaced  with  rigid  reflecting  container  walls. 

Although,  for  testing  purposes,  the  boundary  conditions  in  the  direction  of  uniaxial  flow 
(Z)  remain  as  in  MNCELL  above,  we  note  at  this  stage  that,  one  of  the  rigid  boundaries  in  the 
direction  perpendicular  to  the  slah  can  he  removed  such  that  in  the  presence  of  a  second 
perpendicular  external  force  (e.g  gravity)  or  an  adhesive.  Under  these  conditions  open  surfaces 
can  be  used  for  thin  film  production.  In  this  way,  the  program  can  be  easily  modified  to  simulate 
a  real  process  of  laying  down  a  continuous  solid  tape,  under  gravity,  on  an  adhesive  substrate. 

The  results  for  a  set  of  test  runs  for  the  program  MNSLAB  with  n=12  and  m=6,  and 
density  and  kinetic  energy  (“temperature”)  set  equal  to  1  as  above,  are  summarized  in  Table  3 
below. 


TABLE  3 :  Test  runs  for  MNSLAB  rho=Nro^A^=  1 .0/  T  =1 .0;  slab/rigid  Y 


Run/OP  file 

N 

Time/dt 

Fd 

Vb 

<U> 

<p> 

Comments 

MNOPSLl 

1024 

100000 

10.0 

0.05 

-3.169 

1.574 

1x0.125x4  Equilibration 

MNOPSL2 

1024 

1000000 

10.0 

0.02 

-3.235 

53.12 

2-phase/order  solid 

MNOPSL3 

1024 

1000000 

10.0 

0.5 

-3.415 

1.269 

1 -phase  amorphous 

MNOPSL4 

1024 

100000 

10.0 

0.  1 

-3.174 

1.529 

1 -phase  amorphous 

MNOPSL5 

1024 

100000 

10.0 

0.02 

2.339 

45.33 

2-phase  amorphous  sol. 

MNOPSL6 

1024 

1000000 

10.0 

0.02 

3.392 

54.24 

vac/cryst  zero  flow 

MNOPSL7 

1024 

1000000 

10.0 

0.01 

9.730 

79.62 

Eq.  1x0.625x8  2-phase  cryst 

MNOPSL8 

1024 

1000000 

1.0 

0.01 

-2.429 

12.06 

“  vac/cryst(  1.64)  zero  flow 

MNOPSL9 

1024 

1000000 

1.0 

0.02 

-2.393 

12.32 

“  vac/cryst(l. 64)  zero  flow 
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MNOPSLIO 

1024 

1000000 

0.5 

0.2 

-2.849 

6.602 

“  gas(0.012)/cryst(1.64) 

MNOPSLll 

1024 

1000000 

1.0 

0.1 

-2.427 

11.35 

“  gas(0.001)/cryst(1.64) 

MNOPSL12 

1024 

1000000 

0.5 

0.4 

-2.939 

6.098 

“  gas(~0.065)/  cryst(0.33) 

The  results  show  that  the  same  pattern  of  2-phase  behavior  can  be  obtained  for  the  foil 
periodic  system.  At  very  low  velocities  the  system  crystallizes  completely  such  that  there  is  no 
perceptible  feed  phase,  i.e.  it  appears  t  be  a  vacuum.  There  appears  to  be  a  fairy  narrow  range  of 
2-phase  behavior  ,  when  Fd=10,  because  for  Vb  >0.1,  only  a  single  amorphous  phase  is  obtained. 
It  is  clear  that  the  phase  behavior  depends  strongly  on  the  dimensions  of  the  slab.  For  a  different 
shaped  slab,  which  is  only  half  as  thin  (Runs  6-12),  the  results  are  quite  different. 

c)  lONCELL.f 

A  third  modification  program  has  been  developed  which  can  be  used  for  binary  systems 
and  ionic  systems,  in  this  case  the  pair  potential  takes  the  form  (m=l) 

Eo  ~  Enon-bonded  particleas  ^AB  {  (foAB  A)  — 2ZiZj(  roAB  A)  } 

the  program  has  been  modified  and  placed  on  Orion  ready  to  be  tested  both  for  foil  periodicity, 
and  thin  film  production,  for  an  ionic  salt  such  as  the  solidification  of  sodium  chloride. 

4)  Colloidal  rheometry 

Computational  fluid  mechanics  requires  just  two  ingredients  (i)  solution  of  the 
conservation  equations  of  continuum  motion,  and  (ii)  constitutive  rheological  relations  for  the 
material  constituent  fluids  properties.  The  latter  are  unknown  and  inaccessible  for  complex 
colloidal  fluids,  such  as  powders,  suspensions  of  powders,  plastisols  and  bio-fluids  e.g.  blood. 
Conventional  rheometry,  that  senses  stress  or  strain  rate  at  driving  surfaces,  cannot  detect  the 
spatial  variations  in  velocity  gradients  and  particle  concentration  caused  by  the  shear  field. 

As  demonstrated  already  in  these  computations,  in  principle,  the  obstacles  of 
conventional  rheometry  can  be  eliminated  with  the  use  of  surface-free  periodic  boundaries  and 
the  imposition  of  linear  profiles  artificial  “thermostatting”.  The  situation  at  present,  however,  is 
that  we  are  still  someway  from  obtaining  the  overall  constitutive  rheological  relations  that  will 
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enable  flow  prediction  of  real  colloidal  suspension.  Also,  progress  is  hampered  due  to  a  lack  of 
requisite  quantitative  information  of  the  interparticulate  forces. 

Very  recently,  new  non-intrusive  experimental  methods  are  being  researched  and 
developed,  that  can  sense  the  velocity  gradients,  and  flow  profiles,  and  hence  access  the  stresses, 
and  the  constitutive  rheological  properties  of  complex  non-Newtonian  fluids  in  flow.  Such 
techniques  as  laser  Doppler  anemometry,  NMR  tomography,  positron  tracking  (for  doped 
systems).  All  these  methods  have  inherent  disadvantages.  A  new,  promising  technique  at  an 
early  stage  of  development  which  is  non-intrusive,  relatively  inexpensive  and  fast,  uses  Raman 
light  scattering. 

The  problems  with  conventional  rheometry  for  suspensions  can  be  summarized  as 
follows.  The  thermodynamic  state  of  the  suspension  is  itself  a  function  of  the  rate  of 
deformation;  it  is  this  aspect,  the  lack  of  knowledge  or  information  of  the  change  in  osmotic 
state  and  particle  transport  properties  with  deformation  rate,  that  renders  conventional  rheometry 
inapplicable. 

Consider  a  particulate  fluid  in  steady  shear  flow  in  a  typical  rheometer.  When  steady-state  is 
reached  there  are  forces  at  play  that  make  the  particulate  system  inhomogeneous.  It  follows  that 
the  rheological  information  that  is  obtained  is  peculiar  to  the  geometry  of  the  flow  device.  The 
three  most  common  such  devices  are  the  coaxial  cylinder,  the  cone-and  plate  and  the  capillary 
rheometer. 

When  a  colloidal  suspension  is  at  rest  the  particles  can  come  to  equilibrium  via  Brownian 
motion,  and  they  exhibit  osmotic  properties  of  classical  thermodynamics.  Osmotic  pressures  of 
suspensions  of  colloidal  particles,  of-the-order  1  micron  or  larger,  are  very  small  [ca  lO"^  atm.  at 
f  =  50%  solids].  When  a  system  of  colloidal  micro-particles,  however,  is  subjected  to  shear  flow, 
the  hydrodynamic  forces  on  the  particle  cause  drastic  increases  in  osmotic  pressure  with  the 
square  of  shear  gradient.  At  extreme  rates  of  shear,  for  example  in  very  fast  roll  coating 
operations,  the  osmotic  pressure  can  exceed  one  atmosphere  causing  an  apparent  surface  drying 
effect,  which  is  a  well  known  phenomenon  for  “dilatant”  dense  suspensions,  and  has  to  be 
avoided  by  working  at  shear  rates  lower  than  the  critical  shear  rate  for  dilatancy. 

As  a  consequence  of  the  above  effects,  dense  suspensions  can  become  inhomogeneous 
during  flow  and  processing;  the  particles  tend  to  migrate  fi-om  regions  of  high  shear  to  regions  of 
low  shear.  The  net  result  of  this  readjustment  in  the  three  rheometers  above  is  spurious 
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rheological  data,  because  of  the  incorrect  assumption  of  homogeneity.  In  the  co-axial  cylinder 
the  particles  move  away  from  the  driving  surface  and  towards  the  immobile  surface  (the  cup).  In 
the  cone  and  plate  the  particles  can  cake  on  the  plate,  and  the  effect  in  the  case  of  the  capillary 
rheometer  is  that  particle  move  from  the  pipe  wall,  where  the  shear  rate  is  greatest  towards  the 
pipe  axis  where  the  shear  rate  is  zero.  In  an  extreme  case  of  the  latter  “plug”  flow  results.  This 
effect  in  a  capillary  rheometer  would  cause  the  apparent  shear  viscosity  to  appear  to  be  much 
lower  than  is  would  be  for  the  homogeneous  suspension. 

All  these  undesirable  effects  can  also  occur  in  transport  and  processing  devices  and  need 
to  be  predicted  and  controlled. 

Not  only  is  the  composition  of  the  suspension  in  shear  flow  rheometers  inhomogeneous, 
but  so  also  in  the  shear  gradient.  Under  inhomogeneous  conditions  the  system  further  adjusts  to 
a  high  velocity  gradient  where  the  viscosity  is  low,  in  regions  of  low  concentration,  and  a  low 
gradient,  or  flow  rate,  in  regions  of  high  concentration.  The  extreme  consequence  of  the  above 
inhomogeneities  in  concentration  and  shear  gradient  is  a  problem  well-known  to  chemical 
engineers,  accumulation  and  “caking”  of  the  particles  in  a  static  region  of  zero  flow. 

It  follows  that  the  constitutive  rheological  relationships  for  the  effect  of  a  shear  flow  field 
on  the  stresses,  not  just  the  non-Newtonian  viscosity,  but  also  osmotic  pressures  and  particle 
diffusivities  of  colloidal  suspensions,  are  inaccessible  by  conventional  rheometric  devices. 

On  the  experimental  side,  the  use  of  conventional  rheometry  for  particulate  systems  is 
now  accepted  to  be  untrustworthy.  Real  progress  towards  accurate,  though  perhaps,  more 
limited,  rheological  properties,  is  now  being  made  with  the  development  of  powerful  non- 
intrusive  experimental  techniques  for  obtaining  rheological  information  from  suspension  flow. 
All  have  inherent  advantages  and  disadvantages  regarding  cost,  sample  requirements,  and 
accessibility.  Eventually,  it  is  hoped  that  a  relatively  simple  cost-effective  technique  will 
emerge,  perhaps  the  Ramascope  for  particles  that  it  can  “see”  is  the  most  promising  but  it  is  at 
the  earliest  stage  of  development. 

It  is  to  facilitate  the  development  of  this  non-intrusive  technique  that  a  standard 
crystallizable  monodisperse  colloid  that  is  reproducible  is  required.  To  this  end,  we  have 
explored  a  recent  experimental  prescription  [8,9]  that  promised  to  synthesize,  to  size  order, 
samples  of  monodisperse  latex  with  neutral  surfaces  (hard-spheres)  that  suspend  in  water  with 
neutral  buoyancy. 
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5)  Monodisperse  submicron  latex:  experimental  synthesis 


Starting  from  the  work  in  1982  of  Almog  and  co-workers  [7]  on  the  production  of 
spheres  by  a  single-step  dispersion  polymerization  process.  The  simplicity  of  this  process  was  a 
factor  in  its  popularity  for  study  and  variations  were  made  including  varying  the  monomer, 
studying  the  roughness  of  the  particles  where  the  properties  of  the  steric  stabilizer  were 
investigated,  and  investigating  the  deformabiltiy  of  particles  where  cross  linked  spheres  were 
produced  without  stabilizers.  The  mechanism  of  dispersion  polymerization  is  more  complex  and 
remains  less  well  understood  than  emulsion  polymerization.  As  a  consequence  of  this  ignorance 
It  has  been  difficult,  imtil  recently,  to  prepare  submicron  monodisperse  latex,  to  a  prescribed 
particle  diameter  and  sharp  monodispersity. 

To  produce  particles  in  the  sub-micron  range,  micro-emulsion  polymerization  is  used.  A 
two-stage  process,  micro-emulsion  preparation  and  then  polymerization  usually  achieve  this.  In 
this  system  it  is  accepted  that  the  principle  loci  of  polymerization  are  the  monomer  swollen  latex 
particles.  However  the  particles  produced  by  this  method  exhibited  a  wide  size  distribution 
skewed  to  the  smaller  size  range.  Thus  if  a  narrow  size  range  is  desired  it  is  necessary  to  shorten 
the  nucleation  period  and  avoid  homogeneous  nucleation  in  the  aqueous  phase.  Using  a  redox 
initiator  system,  with  a  faster  dissociation  rate  than  a  thermal  system  lowers  the  nucleation 
period  [references  papers79,80]  this  has  been  coupled  with  the  use  of  a  saw-toothed  blade  mixer 
for  improved  agitation  .  Particles  produced  by  this  method  have  a  very  narrow  distribution  and  a 
size  range  of  100-150  nm.  We  require  particles  around  the  500nm  range. 

Conventional  surfactants  are  small  and  mobile  and  can  migrate  to  the  surface  layer  of  a 
polymeric  film.  This  can  have  a  negative  effect  on  the  application  properties.  One  solution  is  to 
use  a  polymerizable  non-ionic  surfactant  which  can  be  incorporated  into  the  particles  dming 
polymerization  as  described  by  Chem  et  al.  (1996)  [reference  8  paper  81]. 

Turning  now  to  the  very  specific  objective  of  producing  idealized  hard-sphere  particles 
with  the  three  requirements  of  sharp  monodispersity,  a  prescribed  particle  diameter  and  possible 
scale-up  requirement.  Polystyrene  latex  produced  by  dispersion  polymerization  using  the  very 
recent  technique  developed  in  the  last  3  years  in  references  [reference  8  papers  83-85] 
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The  use  of  polyethylene  oxides  by  Lin  and  co-workers  1999  gave  polystyrene  particles  of 
average  size  300  nm.  Liu  and  coworkers  1997  and  1998  have  used  the  same  class  of  stabilizers  in 
dispersion  polymerization  of  styrene  to  give  particles  in  the  range  100  nm  to  1  micron. 
Production  of  particles  “to  order”  at  a  prescribed  diameter  therefore  now  seems  practicable 
subject  only  to  some  preliminary  R&D  work.  The  effect  of  scale-up  is  unpredictable  at  this  stage 
and  needs  to  be  researched. 

The  currently  achievable  standard  deviations  of  the  particle  size  distribution  is  about  1  %. 
Thus  it  should  be  possible,  without  too  much  difficulty,  to  produce  kilogram  quantities  (at  least) 
of  hard-sphere  PS  in  the  size  range  530-535nm  in  the  near  future  if  need  be. 

EXPERIMENT:  METHOD  of  Liu  et  al.  (1998) 

Previously  using  the  Almog  method  monodisperse  polymeric  spheres  in  the  micron  size 
range  were  produced  in  one  step.  These  have  already  been  supplied  to  the  Materials  Research 
Directorate  at  AFRL  [7,8]  Now,  a  new  method  will  be  used  to  produce  monodisperse  polymeric 
spheres  in  the  nanometer  size  range.  This  method  involves  the  dispersion  polymerization  of 
styrene  with  a  small  amount  of  the  polymerizable  w-methoxy  poly-ethylene  oxide  undecyl-a- 
macro-monomer  (PEO-R-MA-40)  as  the  stabilizer. 

We  have  followed  the  prescription  of  Liu  et  al.  (1997).  This  is  based  on  the  dispersion 
polymerization  of  styrene  using  a  small  amount  of  macromonomer  (PEO-R-MA-40)  in  ethanol- 
water  media. 

According  to  Liu  et  al.  the  resultant  size  can  be  formulated.  The  mean  diameter  of  the 
latex  particles  obtained  is  given  by  a  simple  formula  that  depends  on  the  relative  mass  fractions 
of  stabilizer  (PEO-R-MA-40),  surfactant  (AIBN),  and  styrene,  at  the  outset  of  the  polymerization 
as  follows. 

D  X  CF^  =  [PEO-R-MA-40]  ^  [styrene]  ’  [ABN] 

Where  D  is  the  requisite  mean  particle  diameter  and  CF  is  the  fractional  conversion  of  the 
monomer  styrene  to  polymer.  In  the  first  reproducibility  test  of  the  above  the  [styrene]  and 
surfactant  [AEBN]  are  fixed  and  [[PEO-R-MA-40]  varied  such  that  when  the  polymerization  is 
stopped  at  a  given  conversion,  particles  of  size  less  than  SOOnm  are  obtained. 
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The  first  phase  of  this  work  is  completed,  and  has  been  concerned  with  the  synthesis  of 
the  macromonomer  (PEO-R-MA-40),  this  has  successfully  been  done  following  the  prescription 
given  by  Liu  et  al.  This  is  being  used  together  with  styrene  in  the  second  phase  of  the  work  to 
produce  monodisperse  spheres.  The  method  of  Liu  et  al  has  now  been  successfully  reproduced  to 
produce  particles  in  that  size  range.  The  first  product  showed  sub  500nm  particles  with  rather  a 
wider  polydispersity  than  expected,  but  the  samples  have  not  yet  been  fully  analyzed  for  size 
distribution  except  under  a  microscope.  This  research  is  continuing.  A  list  of  all  the  reagents  that 
we  have  had  to  procure  in  order  to  synthesize  500nm  latex  by  this  method  is  summarized  in  the 
appendix  spheres  will  be  available  and  analyzed  before  28/02/01. 

In  summary  we  now  have  the  facility  to  produce  small  laboratory  scale  quantities  of 
monodisperse  aqueous  latex  with  a  mean  particle  size  distribution  of  the  same  range  as  the 
visible  spectrum  of  light.  At  this  stage  it  is  not  possible  to  say  just  how  precisely  a  particular 
color  e.g  green  532  +5  ,  say,  can  be  targeted. 

To  date  we  have  found  no  literature  or  published  attempts  to  scale  up  beyond  the 
laboratory  scale  of  less  than  Order  1kg,  i.e.  to  produce  commercial  quantities.  The  research 
papers  of  Liu  et  al.,  as  with  previous  papers  on  laboratory  dispersion  polymerization,  give  no 
information  on  the  effect  of  scale  up  on  their  formula  that  holds  for  sub-kilogram  quantities  on 
the  1  litre  laboratory  batch  scale. 
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APPENDIX  1 


MNSLAB.f  annotated 


PROGRAM  MNSLAB 
C  MNSLAB  LAST  REVISION  17  NOV  00 

DIMENSION  RX(12000),  RY(12000),  RZ(12000)  particle  positions 

DIMENSION  VX(12000),  VY(12000),  VZ(12000)  particle  velocities 

DIMENSION  FX(12000),  FY(12000),  FZ(12000)  particle  forces 

DIMENSION  PZZ(12000),PS(12000),  SZ(12000)  particle  pressures  and  stresses 

DIMENSION  DISPX(12000),DISPY(12000),DISPZ(12000)  particle  displacements 

DIMENSION  X0(6),  Y0(6),  Z0(6)  initial  lattice  coordinates  (4  for  fee  cell) 

DIMENSION  NIN(8000X  L0(8000),  LI  (8000)  no.  if  particles,  first  and  last  in  n-th  cell 

DIMENSION  DEN(200),VLZ(200),TEE(200),STR(200),POT(200),ZZP(200)  Z-profiles  of  density,  velocity, kinetic 

energy,  stress,  potential  energy  and 


pressure  pzz 

DIMENSION  DXSQ(12),DYSQ(12),DZSQ(12)  mean-square  displacement  by  layer 

DIMENSION  TEELAY(12),POTLAY(12),STRLAY(12)  KE,  PE  and  stress  by  layer 

DIMENSION  XNORM(12),XSUM(12)  accumulators  of  layer  number  densities 


OPEN(5  ,FILE='MND  ATA')  input  data  parameters 

OPEN(6,FILE=’MNRES')  print  out  results 

OPEN(7,FILE='MN  1 024')  configuration  of  positions  and  velocities/start/end 


C 

C  PAIR  POTENTIAL  PHI  =  Ar**(-m)  +  Br**(-n) 

C  PARAMETERS 

READ(5 , 1 00)N,M,CO  attractive  and  repulsive  exponents,  cut-off  in  units  of  ro 

100  FORMAT(2I3,F5.2) 

A  =  FLOAT(N)/FLOAT(M  -  N) 

B  =  FLOAT(M)/FLOAT(N  -  M) 

WRITE(6,200)  N,  M,  A,  B,  CO 
M2  =  M/2 
N2  =  N/2 


C 

C  SYSTEM  STATE 

RE AD(5 ,101  )RHO,TI  system  density  and  initial  temperature 

101  FORMAT(2F8.4) 

C  STEADY  STATE  PROCESS  VARIABLES 

READ(5,102)FDZ,VBZ,NVB  system  driving  force,  boundary  velocity  (compaction  velocity),  NVB  is  z-boundary 
condiotion 

102  FORMAT(2F8.4,I3) 

IF(NVB.EQ.  1  )WRITE(6,22 1 ) 

IF(NVB.EQ.2)WRITE(6,222) 

IF(NVB.EQ.3)WRITE(6,223) 

C 

WRITE  (6,201)  RHO,  TI,  FDZ,  VBZ 
C 

C  PROGRAM  CONSTANTS 

READ(5 ,103)  NA  total  number  of  particles 

103  FORMAT(I8) 

WRITE  (6,202)  NA 
XNA=NA 

READ(5,104)NRUN  index  of  run  0,l,2,3,etc.  lattice  start  for  runO 

104  FORMAT(I3) 

READ(5,105)  NDT,  DT,  ISCALE  number  of  time  steps,  time  increment,  thermostat  or  not  (0  or  1) 

105  FORMAT(  18,  F8.4, 13) 

C  PRINT  FREQUENCY 

READ(5,106)NSUB 

106  FORMAT(I8) 

XNSUB  =  NSUB 

WRITE  (6,203)  NRUN,  NDT,  NSUB,  DT,  ISCALE 
C 

C  BOX  SHAPE  DIMENSIONS  volume  in  reduced  units  XF  x  YF  x  ZF 
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XF=  1.0 

READ(5,107)YF,ZF 

107  FORMAT(2F8.4) 

WRITE(6,206)XF,YF,ZF 

C 

C  CONVERSIONS  TO  PROGRAM  UNITS 

C  LENGTH  is  length  of  box  in  X,  MASS  is  particle  mass,  and  TIME  is  DT 
XL=(XNA/(RHO*YF*ZF))**(l  .0/3.0) 

RO  =  XF/XL  ^stance  of  minimum  potential  in  program  units 

EPS=R0*R0*DT*DT  minimum  pair  energy  in  program  units 
A=A*EPS*(R0**M) 

B=B*EPS*(R0**N) 

AM=A*M 
BN=B*N 
CO=CO*RO 
C02-C0*C0 
XDIM  =  XF/RO 
YDIM  =  YF/RO 
ZDIM=  ZF/RO 

WRITE(6,213)  XDIM,  YDIM,  ZDIM 
FDZ-FDZ*DT*DT*R0 
VBZ  =  VBZ*DT*R0 
C 

C  PROFILE  HISTOGRAM  PARAMETERS 
READ(5, 1 08)NLAY,XHIST 

108  FORMAT(I3,F6.1) 

XNLAY  =  NLAY 
SPAC=ZF/XNLAY 

VOLHIS  =  (ZF/XHIST)  *  YF  *  XF/(R0**3) 

NHIST  -  IFIX(XHIST+0. 1 ) 

C 

C  LINK  CELLS 

READ(5, 1 09)  NLX,  NLY,  NLZ  numbers  of  cells  in  each  direction 

109  FORMAT(3I3) 

WRITE(6,204)  NLX,  NLY,  NLZ 
NL2-NLX*NLY 
NL3-NLX*NLY*NLZ 
XNL  =  NLX 

YNL  =  NLY 
ZNL  =  NLZ 

CMAXX  =  XF/(XNL  *  RO)  maximum  spherical  pair  potential  cut-off  allowed 
CMAXY  =  YF/(YNL  *  RO) 

CMAXZ  =  ZF/(ZNL*R0) 

WRITE  (6,2 1 0)  CMAXX,CMAXY,CMAXZ  cot  off  must  be  less  than  the  smallest  of  these 
C 

IF(NRUN.GE.l)GO  TO  2 

C  LATTICE  CONFIGURATION  AND  INITIAL  VELOCITIES 
READ  (5,1 10)  NCA,NXA,NYA,NZA 

110  FORMAT(4I3) 

WRITE(6,209)  NCA,NXA,NYA,NZA 
XNXA  =  NXA/XF 
XNYA  =  NYA/YF 
XNZA  =  NZA/ZF 
1=1 

TEMPX  =  0.0 
TEMPY  =  0.0 
TEMPZ  =  0.0 
D0  3  J=1,NCA 

READ(5,1 1 1)  X0(J),Y0(J),Z0(J) 

111  FORMAT(3F8.4) 

WRITE(6,21 1)  J,  X0(J),Y0(J),Z0(J) 

D0  3IX=1,NXA 
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D0  3IY=1,NYA 
D0  3IZ=1,NZA 

C  initial  positions  for  a  face  centered  cubic  lattice 
RX(I)=((IX- 1  )+X0(  J))/XNXA 
RY(I)=((IY- 1  )4-Y0(  J))/XNYA 
RZ(I)=((IZ- 1  )+Z0(  J))/XNZA 

C  use  of  initial  positions  to  generate  initial  random  velocities  (caution-equilibrate  well) 

VY(I)  =  RX(I)  -  0.5*XF 
VZ(I)  =  RY(I)  -  0.5  *XF 
VX(I)  =  RZ(I)  -  0.5*ZF 
C  initial  temperatures/kinetic  energies 
TEMPX=TEMPX  +  VX(I)**2 
TEMPY^TEMPY  +  VY(I)**2 
TEMPZ=TEMPZ  +  VZ(I)**2 
1=1+1 

3  CONTINUE 

TEMPX  =  TEMPX/(XNA  *  EPS) 

TEMPY  =  TEMPY/(XNA  *  EPS) 

TEMPZ  =  TEMPZ/(XNA  *  EPS) 

C  scaling  factor  required  to  renormalise  temp  in  each  direction  separately 
FACX  =  SQRT(TI/TEMPX) 

FACY  =  SQRT(TI/TEMPY) 

FACZ  =  SQRT(TI/TEMPZ) 

GO  TO  24 
C 

C  READ  IN  CONFIGURATION  FROM  FILE  7 
2  CONTINUE 
FACX  =  1.0 

FACY  =1.0 
FACZ  =1.0 
DO  14 1=1, NA 

READ(7,250)  RX(I),RY(I),RZ(I),VX(I),VY(I),VZ(I) 

C  check  that  Y-  values  (slab  dimension  with  rigid  walls)  are  not  exactly  zero 
IF(RY(I).LT.0.000001)  RY(I)  =  0.000001 
IF(RY(I).GE.YF)  RY(I)  =  YF  -0.000001 
14  CONTINUE 
C 

24  CONTINUE 
C  CELL  OCCUPANCY  TABLES 

C  in  this  code  N  particles  are  indexed  according  to  position  starting  with  l-Ll(l)  in  cell(l),  L0(2)-L1(2),  in 

CELL(2)  up  to  L0{N)-n  in  CELL  (NL3) 

J=  1 

DO  20  IP=1  ,NL3  (NL3  is  total  number  of  cells  NLX*NLY*NLZ) 

LO(IP)  =  J 
L1(IP)  =  0 
DO  21  1=  1,NA 

IX  =  INT(  RX(I)  *  XNL/XF  +  1 .0) 
lY  =  INT(  RY(I)  *  YNLA^F) 

IZ  =  INT(RZ(I)*ZNL/ZF) 

ICI  =  DC  +  NLX*IY  +  NL2*IZ 
IF(ICLNE.IP)  GO  TO  21 
FX(J)=RX(I) 

FY(J)=RY(I) 

FZ(J)=RZ(I) 

L1(IP)  =  J 
J  =  J+  1 

21  CONTINUE 

IF(L1(IP).EQ.0)L0(IP)  =  0 
20  CONTINUE 
C 

DO  96  1=1, NA 
RX(I)=FX(I) 
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RY(I)=FY(I) 

RZ(I)=FZ(I) 

96  CONTINUE 
C 

C  WRITE  OUT  CELL  OCCUPANCY 
NCOUNT=0 
3999  CONTINUE 
WRITE(6,208) 

D0  69IZ=1,NLZ 
JZ  =  IZ-1 

DO  69  IY=1,NLY 

JY  =  IY~1 

DO  69  IX  =  1,NLX 

P=  DC  +  NLX  *  JY  +  NL2  *  JZ 

NIN(JP)=L1  (JP)-L0(P)+ 1 

IF(L0(JP).EQ.0)NIN(P)  =  0 

WRITE(6,*)JP,L0(JP),L1(JP),NIN(JP) 

69  CONTINUE 
C 

C  CHECK  CELL  OCCUPANCY  INDEX 
D0  31J=1,NA 

IX  =  INT(  RX(J)  *  XNL/XF  +  1 .0) 
lY  =  INT(  RY(J)  *  YNLA^F) 

IZ  =  INT(  RZ(J)  *  ZNL/ZF) 

JP  -  IX  +  NLX*IY  +  NL2*IZ 
IF(J.GE.LO(JP).ANDJ.LE.L1(P))GO  TO  31 
WRITE  (6,220)  J,JP,L0(P),L1(JP) 

31  CONTINUE 
IF(NCOUNT.NE.O)  GO  TO  998 

C 

C  WRITE  INITIAL  CONFIGURATION 
TEMPX  =  0.0 
TEMPY  =  0.0 
TEMPZ  =  0.0 
WRITE  (6,214) 

D0  32W,NA 

IF(LLE.  1 0)WRITE(6,250)  RX(I),RY(I),RZ(I),VX(I),VY(I),VZ(I) 

VX(I)  =  VX(I)*FACX 

VY(I)  =  VY(I)*FACY 

VZ(I)  =  VZ(I)*FACZ 

TEMPX=TEMPX  +  VX(I)**2 

TEMPY=TEMPY  +  VY(I)**2 

TEMPZ^TEMPZ  +  VZ(I)**2 

32  CONTINUE 

TEMPX  =  TEMPX/(XNA  *  EPS) 

TEMPY  =  TEMPY/pCNA  *  EPS) 

TEMPZ  =  TEMPZ/(XNA  *  EPS) 

TEMP  =  (TEMPX  +  TEMPY  +  TEMPZ)/3.0 
WRITE(6,207)  TEMP,  TEMPX,  TEMPY,  TEMPZ 
C 

C  INITIAL  MOMENTA 
PX=0.0 
PY=0.0 
PZ=0.0 
DO  5  1=1, NA 
PX=PX+VX(I) 

PY=PY+VY(I) 

PZ=PZ+VZ(I) 

5  CONTINUE 
C 

C  ZERO  MOMENTA 
X=PX/XNA 
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CJU  O  cn<nUO 


Y=PY/XNA 

Z=PZ/XNA 

TEMPX  =  0.0 
TEMPY  =  0.0 
TEMPZ  =  0.0 
DO  6  1=1, NA 
VX(I)=VX(I)-X 
VY(I)=VY(I)-Y 
VZ(I)=VZ(I)-Z 

TEMPX=TEMPX  +  VX(I)**2 
TEMPY=TEMPY  +  VY(I)**2 
TEMPZ=TEMPZ  +  VZ(I)**2 
6  CONTINUE 

TEMPX  =  TEMPX/pCNA  *  EPS) 

TEMPY  =  TEMPY/(XNA  *  EPS) 

TEMPZ  =  TEMPZ/(XNA  *  EPS) 

ZERO  ALL  ACCUMULATORS 
NCOUNT  =  0 
NPRINT  =  0 
ZNFLOW  =  0.0 
TEMPAV  =  0.0 
TEMXAV  =  0.0 
TEMYAV  =  0.0 
TEMZAV  =  0.0 
PVAV  =0.0 
UAV  =0.0 
PHIAV  =  0.0 
PSIAV  =  0.0 
PVSQ  =0.0 
PHISQ  =  0.0 
PSISQ  =  0.0 
PSIPHI  =  0.0 
SET  ARRAYS  TO  ZERO 
D0  7L=  1,  100 
DEN(L)  =  0.0 

VLZ(L)  =  0.0 
TEE(L)  =  0.0 
STR(L)  =  0.0 
POT(L)  =  0.0 
ZZP(L)  =  0.0 
CONTINUE 
DO  25L=  1,NLAY 
XNORM(L)  =  0.0 
TEELAY(L)  =  0.0 
POTLAY(L)  =  0.0 
STRLAY(L)  =  0.0 
;  CONTINUE 
DO  261=  1,NA 

DISPX(I)  =  0.0 
DISPY(I)  =  0.0 
DISPZ(I)  =  0.0 
)  CONTINUE 
ZNFLOW  =  0.0 

ZNFLOW  =  0.0 

SUM  FORCE  ON  EACH  ATOM 
12  CONTINUE 
C 

C  THE  VELOCITIES  ARE  SCALED  IF  ISCALE=1 
C  EACH  DIRECTION  SEPERATELY 
FACX  =  SQRT(TI/TEMPX) 
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FACY  =  SQRT(TI/TEMPY) 
FACZ  =  SQRT(TI/TEMPZ) 

IF  (ISCALE.EQ.  1)  GOTO  46 
IF  (NCOUNT.EQ.O)  GOTO  46 
FACX=1.0 
FACY  =1.0 
FACZ  =1.0 
C 

46  CONTINUE 
DO  8  1=1, NA 
PZZ(I)  =  0.0 
SZ(I)  =  0.0 
PS(I)  =  0.0 

VX(I)  =  VX(I)*FACX 
VY(I)  =  VY(I)*FACY 
VZ(I)  =  VZ(I)*FACZ 
8  CONTINUE 
PHIM=0.0 
PHIN=0.0 
IP=0 

D0  9IZ=1,NLZ 
D0  9IY=1,NLY 
D0  9IX=1,NLX 


IP  =  IP  +  1 

C  . CELL  IP  EMPTY 

IF  (LO(IP).EQ.O)  GOTO  39 
C  . CALCULATE  JP 


JY=1 

DO  40  KX  =  IX-1,IX+1 
DO  40  KZ  =  IZ,IZ+1 

IF(KX.EQ.IX-l.AND.KZ.EQ.IZ)GO  TO  40 
CX  =  0.0 
CY  =  0.0 
CZ  =  0.0 
JX  =  KX 
JZ  =  KZ 

C  . PERIODIC  BOUNDARY  CONDITIONS 

C  this  is  a  slab  with  rigid  walls  in  the  Y-direction  and  periodicity  in  X  and  Z  only 
IF(JX-NLX4)  73,70,73 
70  JX=1 

CX  =  XF 

73  IF(JZ-NLZ-l)  75,74,75 

74  JZ  =  1 
CZ  =  ZF 

75  IF  (JX)  79,  76,  79 

76  JX  =  NLX 
CX  =  -XF 

79  IF  (JZ)  82,  80,  82 

80  JZ  =  NLZ 
CZ  =  -ZF 

C 

82  JP  =  JX  +  (JY-1)*NLX  +  (JZ-1)*NL2 
IF(L0(JP).EQ.0)  GO  TO  41 
C 

N1=L0(IP) 

N2=L1(IP) 

D019I=N1,N2 

C 

XR=RX(I) 

YR=RY(I) 

ZR=RZ(I) 

C 
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N3=L0(JP) 

N4=L1(JP) 

DO  18  J=N3,N4 

IF(IP.EQ.JP.AND.LGE.J)GO  TO  18 
C  calculate  forces  between  particles  I  (I)  and  j(J) 
X=RX(J)>XR+CX 
Y=RY(J)-YR+CY 
Z=RZ(J)-ZR+CZ 
XX-X*X 
YY  =  Y*Y 
ZZ  =  Z*Z 

R2-XX  +  YY  +  ZZ 
IF(R2.GT.C02)  GOTO  18 
R2=L0D0/R2 
RM=R2*R2*R2 
RN=RM*RM 
PHIN  =  PHIN  +  RN 
PHIM  -  PHIM  +  RM 
F  =  A*RM  +  B*RN 
PS(I)-PS(I)  +  F 
PS(J)  =  PS(J)  +  F 
C 

F  =  AM*RM  +  BN*RN 

F=F*R2 

X=F*X 

Y=F*Y 

Z=F*Z 

FZZ  =  ZZ  *  F 
SZZ  =  PCX-ZZ)  *  F 
VX(I)  =  VX(I)  -  X 
VY(I)  =  VY(I)  -  Y 
VZ(I)-VZ(I)-Z 
VX(J)  =  VX(J)  +  X 
VY(J)  =  VY(J)  +  Y 
VZ(J)  =  VZ(J)  +  Z 
PZZ(I)-PZZ(I)-FZZ 
PZZ(J)  =  PZZ(J)  -  FZZ 
SZ(I)  =  SZ(I)  +  SZZ 
SZ(J)  =  SZ(J)  +  SZZ 

18  CONTINUE 

19  CONTINUE 
41  CONTINUE 
40  CONTINUE 
39  CONTINUE 
9  CONTINUE 
C 

NCOUNT  =  NCOUNT+  1 
XCOUNT  =  NCOUNT 
C  UPDATE  CONFIGURATION 
C 

TEMP  =  0.0 
TEMPX=  0.0 
TEMPY^O.O 
TEMPZ=  0.0 
PX=0.0 
PY=0.0 
PZ=0.0 

DO  27L=  1,NLAY 

DXSQ(L)=0.0 

DYSQ(L)=0.0 

DZSQ(L)=0.0 

XSUM(L)-0.0 
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27  CONTINUE 
C 

DO  10I=1,NA 
C  SUM  MOMENTA 
PX=PX+VX(I) 

PY=PY+VY(I) 

PZ=PZ+VZ(I) 

C  ACCUMULATE  Z-PROFILES 
K  =  INT(RZ(I)  *  XHIST/ZF)  +  1 
DEN(K)  =  DEN(K)+  1.0 

VLZ(K)  =  VLZ(K)  +  VZ(I) 

STR(K)  =  STR(K)  +  SZ(I) 

ZZP(K)  =  ZZP(K)  +  PZZ(I) 

POT(K)  =  POT(K)  +  PS(I) 

C  LAYER  PROPERTIES 

L  =  INT(RZ(I)  *  XNLAY/ZF)  +  1 
C  DIFFUSION 

DISPX(I)  =  DISPX(I)  +  VX(I) 

DISPY(I)  =  DISPY(I)  +  VY(I) 

DISPZ(I)  =  DISPZ(I)  +  VZ(I) 

DXSQ(L)  =  DXSQ(L)  +  DISPX(I)*DISPX(I) 

DYSQ(L)  =  DYSQ(L)  +  DISPY(I)*DISPY(I) 
DZSQ(L)  =  DZSQ(L)  +  DISPZ(I)*DISPZ(I) 
XSUM(L)  =  XSUM(L)  +  1.0 


C 

EKIX  =  VX(I)**2 
EKIY  =  VY(I)**2 
EKIZ  =  VZ(I)**2 

EKI  =  (EKIX  +  EKIY  +  EKIZ)/(3.0  *  EPS) 
TEE(K)  =  TEE(K)  +  EKI 
TEMPX  =  TEMPX  +  EKIX/(XNA  *  EPS) 
TEMPY  =  TEMPY  +  EKIY/(XNA  *  EPS) 
TEMPZ  =  TEMPZ  +  EKIZ/pCNA  *  EPS) 
TEMP  =  TEMP  +  EKI/XNA 
TEELAY(L)  =  TEELAY(L)  +  EKI 
POTLAY(L)  =  POTLAY(L)  +  PS(I) 
STRLAY(L)  =  STRLAY(L)  +  SZ(I) 
XNORM(L)  =  XNORM(L)  +  1 .0 
C  OLD  CELL 

IX  =  INT(  RX(I)  ♦  XNL/XF  +  1.0) 
lY  =  INT(  RY(I)  *  YNLAT) 

IZ  =  INT(  RZ(I)  *  ZNL/ZF) 

ICI  =  IX  +  NLX*IY  +  NL2*IZ 
C  EXTERNAL  FORCE 

VZ(I)  =  VZ(I)-FDZ 
C  NEW  POSITIONS 
RX(I)  =  RX(I)  +  VX(I) 

RY(I)  =  RY(I)  +  VY(I) 

RZ(I)  =  RZ(I)  +  VZ(I) 

C  PERIODIC  BOUNDARY  IN  X 

IF(RX(I).GE.XF)  RX(I)  =  RX(I)-XF 
IF(RX(I).LT.0.0)  RX(I)  =  RX(I)  +  XF 
C  RIGID  BOUNDARY  IN  Y 

IF(RY(I).GE.YF)  RY(I)  =  RY(I)  -  VY(I) 
IF(RY(I).GE.YF)  VY(I)  =  -  VY(I) 
IF(RY(I).LT.0.0)  RY(I)  =  RY(I)  -  VY(I) 
IF(RY(I).LT.0.0)  VY(I)  =  -  VY(I) 

C  RIGID  BOUNDARY  AT  Z=ZF 
IF(RZ(I).GE.ZF)  VZ(I)  =  -VZ(I) 
IF(RZ(I).GE.ZF)  RZ(I)  =  ZF  -  0.000001 
IF(RZ(1).GE.0.0)  GO  TO  1000 


C  Z=0  BOUNDARY  CONDITIONS 

CifNVB  l=PERIODIC  2=MOVING  3=RIGID 

IF(NVB.EQ.3)  GO  TO  1001 
IF(NVB.EQ.l)  GO  TO  1002 
C  MOVING  BOUNDARY  AT  Z=0 
IF(-VZ(I).GT.VBZ)  GO  TO  1001 
1002  CONTINUE 
RZ(I)  =  RZ(I)  +  ZF 

C  accumulate  the  flow  rate  in  the  z-direction 
ZNFLOW  =  ZNFLOW  +1.0 
GO  TO  1000 
1001  CONTINUE 
RZ(I)=RZ(I)-VZ(I) 

VZ(I)  =  -VZ(I) 

C  accumulate  the  collision  (reversal)  rate  at  the  +Z  semi-permeablr  boundary. 

ZNCOLL  =  ZNCOLL  +  1 .0 
1000  CONTINUE 
C  NEW  CELL 

IF(RX(I).GE.XF)RX(I)  =  XF  -  0.000001 
IF(RY(I).GE. YF)RY(I)  =  YF  -  0.00000 1 

IF(RZ(I).GE.ZF)RZ(I)  =  ZF  -  0.000001 
IX  =  INT(  RX(I)  *  XNL/XF  +1.0) 
lY  =  INT(  RY(I)  *  YNLA^F) 

IZ  =  INT(  RZ(I)*ZNL/ZF) 

IP  =  IX  +  NLX*IY  +  NL2*IZ 

IF(IP.LE.0.OR.IP.GT.NL3)  GO  TO  999 
C 

C  UPDATE  PARTICLE  INDICES 

C  Once  a  particle  has  moved  across  a  cell  boundary,  the  index  of  all  the  other  particles  changes  to  restore  the  addressing  system 
for  locating  particles  by  cell.  This  code  updates  the  indices  of  all  the  particles. 

IF(IP.EQ.ICI)GO  TO  34 
X=RX(I) 

Y=RY(I) 

Z-RZ(I) 

XX=VX(I) 

YY=VY(I) 

ZZ=VZ(I) 

XXX  =  DISPX(I) 

YYY  -  DISPY(I) 

ZZZ  -  DISPZ(I) 

NIN(ICI)=NIN(ICI)-1 
NIN(IP)  =  NIN(IP)+1 
IF(IP.GT.ICI)GO  TO  1 12 
C 

J  =  L0(ICI) 

IF(L0(ICI).NE.L1  (ICI))GO  TO  1 39 
LO(ICI)  =  0 
Ll(ICI)  =  0 
139  CONTINUE 
L  =  ICI 

DO  117  JP  =  IP+1,ICI 
K  =  L-1 

IF(L0(L).NE.0)  L0(L)=  L0(L)  +  1 

IF(L1(K).NE.0)  L1(K)=  L1(K)  +  1 
IF(L0(L).NE.0)  J=L0(L)-1 
L-L-1 
117  CONTINUE 

IF(L1(IP).NE.0)  GO  TO  127 
LO(IP)  =  J 

L1(IP)  =  J 
127  CONTINUE 
J=I 
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115  K=J-1 

IF(J.LE.L1(IP))G0  TO  33 
RX(J)=RX(K) 

RY(J)-RY(K) 

RZ(J)=RZ(K) 

VX(J)=VX(K) 

VY(J)=VY(K) 

VZ(J)=VZ(K) 

DISPX(J)=DISPX(K) 

DISPY(J)=DISPY(K) 

DISPZ(J)=DISPZ(K) 

J=J-1 

GO  TO  1 15 
C 

112  CONTINUE 
J  =  L1(ICI) 

IF(LO(ICI).NE.L1(ICI))GO  TO  129 
LO(ICI)  =  0 
L1(ICI)  =  0 
129  CONTINUE 
DO  94L  =  ICI,IP-1 
K  =  L+1 

IF(LO(K).NE.O)  L0(K)=  L0(K)  -  1 

IF(L1(L).NE.0)  L1(L)=  L1(L)  -  1 
IF(L1(L).NE.0)  J=L1(L)+1 
94  CONTINUE 

IF(L1(IP).NE.0)  GO  TO  128 
L1(IP)  =  J 
LO(IP)  =  J 
128  CONTINUE 
J=I 

118  K-J+1 

IF(J.GE.LO(IP))GO  TO  33 
RX(J)=RX(K) 

RY(J)=RY(K) 

RZ(J)=RZ(K) 

VX(J)=VX(K) 

VY(J)=VY(K) 

VZ(J)=VZ(K) 

DISPX(J)=DISPX(K) 

DISPY(J)=DISPY(K) 

DISPZ(J)=DISPZ(K) 

J=J+1 

GOTO  118 
C 

33  CONTINUE 
RX(J)=X 
RY(J)=Y 
RZ(J)=Z 
VX(J)=XX 
VY(J)=YY 
VZ(J)=ZZ 

DISPX(J)=XXX 

DISPY(J)=YYY 

DISPZ(J)=ZZZ 

34  CONTINUE 
10  CONTINUE 
C 

C  ACCUMULATE  TIME  AVERAGES 
C 

PHIM=A*PHIM/(EPS*XNA) 

PHIN=B*PHIN/(EPS*XNA) 
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PSI=-(M*PHIM  +  N*PHIN) 

PHI=PHIN+PHIM 
U=  PHI+1.5*TEMP 
PV  =  1.0-PSI/(TEMP*3.0) 

TEMPAV  =  TEMPAV  +  TEMP 
TEMXAV  =  TEMXAV  +  TEMPX 
TEMYAV  =  TEMYAV  +  TEMPY 
TEMZAV  =  TEMZAV  +  TEMPZ 
UAV  =UAV  +U 
PVAV  =PVAV  +PV 
PHIAV  =PHIAV  +PHI 
PSIAV  =  PSIAV  +  PSI 
PHISQ  =  PHISQ  +  PHI  *  PHI 
PSISQ  =  PSISQ  +  PSI  *  PSI 
PSIPHI  =  PSIPHI  +  PSI  *  PHI 
C 

C  REGULAR  PRINT  OUT  EVERY  NPRINT  TIME  STEPS 
C 

IF(NCOUNT.EQ.l)  GO  TO  61 
IF(NSUB.EQ.1)G0T0  61 
IF(NCOUNT.EQ.NPRINT)  GO  TO  61 
GO  TO  60 
61  CONTINUE 

WRITE(6,205)  NCOUNT,  PX,  PY,  PZ,  U,  PHIM,  PHIN,  PV 
WRITE(6,207)  TEMP,  TEMPX,  TEMPY,  TEMPZ 
C 

IF(NRUN.EQ.O)GO  TO  56 

D0  55L=NLAY,1,-1 

XMSD  =  DXSQ(L)/(XSUM(L)*R0*R0) 

YMSD  =  DYSQ(Ly(XSUM(L)*R0*R0) 

ZMSD  =  DZSQ(L)/(XSUM(L)*R0*R0) 

DX  =  XMSD/(2.0*XNSUB*DT) 

DY  =  YMSD/(2.0*XNSUB*DT) 

DZ  =  ZMSD/(2.0*XNSUB*DT) 

WRITE(6,450)L,DX,DY,DZ 

55  CONTINUE 
DO30I-  1,NA 

DISPX(I)  =  0.0 
DISPY(I)  =  0.0 
DISPZ(I)  =  0.0 
30  CONTINUE 
C 

C  RATE  OF  FLUX 

56  CONTINUE 

X  =  ZNFLOW*R0*R0/(XNSUB*DT) 

WRITE(6,2 1 6)ZNFLOW,X,ZNCOLL 
ZNFLOW  =  0.0 
ZNCOLL^O.O 
C 

NPRINT  =  NPRINT  +  NSUB 
C 

60  CONTINUE 

IF(NCOUNT.GE.NDT)  GOTO  1 1 
GOTO  12 
C 

C  END  OF  RUN 
1 1  CONTINUE 
C  FINAL  AVERAGES 

TEMPAV  =  TEMPAV  /  XCOUNT  average  temperature 

TEMXAV  =  TEMXAV  /  XCOUNT 
TEMYAV  =  TEMYAV  /  XCOUNT 
TEMZAV  =  TEMZAV  /  XCOUNT 
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PVAV  =  PVAV/XCOUNT 
UAV  =  UAV/XCOUNT 
PHIAV  =  PHIAV/XCOUNT 
PSIAV  =  PSIAV/XCOUNT 
PVSQ  =  PVSQ/XCOUNT 
PHISQ  =  PHISQ/XCOUNT 
PSISQ  =  PSISQ/XCOUNT 
PSIPHI  -  PSIPHI  /  XCOUNT 
WRITE  (6,355)  NCOUNT 
WRITE  (6,360)  UAV 
WRITE  (6,370)  TEMPAV 
WRITE  (6,380)  PVAV 
WRITE  (6,390)  PHIAV 
WRITE  (6,400)  PSIAV 
WRITE  (6,420)  PHISQ 
WRITE  (6,430)  PSISQ 
WRITE  (6,440)  PSIPHI 
WRITE  (6,460)  TEMXAV 
WRITE  (6,470)  TEMYAV 
WRITE  (6,480)  TEMZAV 
WRITE  (6,280) 

C 

C  NORMALISE  AND  PRINT  OUT  Z-PROFILES 
C 

DO  35L  =  NHIST,1,-1 
Z  =  (ZF/XHIST)  *  (L  -  1)/R0 
D  =  DEN(L)/(XCOUNT  *  VOLHIS) 

V  =  VLZ(L)/(DEN(L)*DT*R0) 

W  =  D*V 

IF  (D.GT.0.00000001)  GOTO  42 
S  =  0.0 
T=0.0 
P  =  0.0 
F  =  0.0 
GOTO  43 

42  S  =  -STR(L)/(XCOUNT  *  2.0  *  EPS  *  VOLHIS) 

T  =  TEE(L)/DEN(L) 

P  =  POT(L)/(  DEN(L)  *  2.0  *  EPS) 

F  =  D  *  T  -  ZZP(L)/(VOLHIS  *  2.0  *  XCOUNT  *  EPS) 

43  WRITE(6,290)L,Z,D,V,W,T,P,S,F 

35  CONTINUE 
C 

C  RMSD  POTENTIAL  AND  KINETIC  ENERGIES  PER  LAYER 
C 

WRITE(6,320) 

D0  36L=NLAY,1,-1 

IF(XSUM(L).EQ.O)  GO  TO  36 
XMSD  =  DXSQ(L)/(XSUM(L)*R0*R0) 

YMSD  -  DYSQ(L)/(XSUM(L)*R0*R0) 

ZMSD  =  DZSQ(L)/(XSUM(L)*R0*R0) 

DX  =  XMSD**0.5 
DY  =  YMSD**0.5 
DZ  =  ZMSD**0.5 

TLAY  =  TEELAY(L)/(XNORM(L)) 

PLAY  =  POTLAY(L)/(XNORM(L)  *  2.0  *  EPS) 

SLAY  =  STRLAY(L)/(XNORM(L)  *  2.0  *  EPS  *  VOLHIS) 
WRITE(6,300)  L,  DX,  DY,  DZ,  TLAY,  PLAY,  SLAY 

36  CONTINUE 

C  WRITE  OUT  FINAL  CONFIGURATION 
REWIND  7 
DO  13  1=1, NA 

IF(RY(I).LT.0.000001)  RY(I)  =  0.000001 


average  pressure  (<pv  /  NkT>) 
average  internal  energy 
average  potential  energy 
average  virial 
mean  square  pressure 
mean  square  potential 
mean  square  virial 
covariance  potential-virial 
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IF(RY(I).GE.YF)  RY(I)  =  YF  -0.000001 
WRITE(7,250)  RX(I),RY(I),RZ(I),VX(I),VY(I),VZ(I) 

13  CONTINUE 
WRITE(6,215) 

DO  37  1=1,50 

WRITE(6,250)  RX(I),RY(I),RZ(I),VX(I),VY(I),VZ(I) 

37  CONTINUE 
C 

C  FORMAT  WRITE  STATEMENTS 
C 

200  FORMAT(/,'n-m  PARAMETERS’, 2(13), 6X, 'A  B’,2(F6.2),2X,’CO=’,F5.2,/) 

201  F0RMAT("RH0=”,F8.4,4X,”TI=”,F8.4,4X;FDZ=’,F8.4,4X,'VBZ=',F8.4,/) 

202  FORMAT(’NUMBER  of  PARTICLES  =',I6,/) 

203  FORMAT(’NRUN=',I3,4X,’NDT=',2I8,4X,’DT=’,F8.4,4X,’ISCALE=’,I3,/) 

204  FORMAT(’NLX=',I2,2X,'NLY=’,I2,2X,'NLZ=’,I2,/) 

205  FORMAT(/,TIME',I8,'  PX  PY  PZ’,3(E9.3),2X,’U',3F7.3,2X,’Z',F8.3) 

206  FORMATCBOX  DIMENSIONS  XF  YF  ZF,3F8.4,/) 

207  FORMAT(’TEMP’,F8.4,2X,’TX  TY  TZ',3(2X,F8.4)) 

208  FORMAT(/, 'CHECK  ON  LINK  CELL  ARRAYS’,/) 

209  FORMAT(/,I3,’  ATOMS  IN  UNIT  CELL’, 313,/) 

210  FORMAT('CMAXXYZare’,3(F8.4)/) 

211  FORMAT('C0-0RDS  XYZ',I2,3F8.4) 

213  FORMAT(  ’XDIM=  •,F8.4,2X,'YDIM=  ’,F8.4,2X,’ZDIM=  ’,F8.4,/) 

214  FORMAT(/, 'INITIAL  CONFIGURATION’,/) 

2 1 5  FORM AT(/,’FIN AL  CONFIGURATION',/) 

216  FORMAT('ZNFLOW=’,F10.0,2X,E1  L4,8X,'ZNCOLL=’,F10.0,/) 

220  FORMATC'OUT  OF  CELL”,I6,3(2X,I5),/) 

22 1  FORM AT(/, 'PERIODIC  BOUNDARY  AT  Z=0',/) 

222  FORMAT(/,'MOVING  BOUNDARY  AT  Z=0',/) 

223  FORMAT(/, 'RIGID  BOUNDARY  AT  Z=0',/) 

250  FORMAT(3F10.6,2X,3E12.4) 

280  FORMAT(/,8X,’Z/ro',4X,'  RHO’,6X,’VZ',6X,’flux’,6X,'  KE’, 

1 6X,'PE',6X,'STRESS',4X,'PZZ’,/) 

290  FORMAT(I4,F8.2,ll(F9.3)) 

300  FORMAT(I5,6(4X,F8.4)) 

320  FORMAT(/, ’LAYER', 6X,’RMSDX’,7X,'RMSDY’,7X,'RMSDZ’,7X,’TEELAY',6X, 
1  ’POTLAY',6X,'STRL  AY’,/) 

355  FORMAT(/,/,’TIME  AVERAGES’, 6X,’NCOUNT=', 18,/) 

360  FORMATC  ENERGY  =  ’,F10.4) 

370  FORMATC  TEMP  =  ',F10.4) 

380  FORMATC  PRESSURE  =  ’,F10.4) 

390  FORMATC  POTENTIAL  =  ',F10.4) 

400  FORMATC  VIRIAL  =  ’,F10.4) 

420  FORMATC  MSQPOT  =  ’,F10.4) 

430  FORMATC  MSQ  VIRIAL  =  ’,F10.4) 

440  FORMATC  Av  PSIPHI  =  ’,F10.4) 

460  FORMATC  TEMP  X  =  ’,F10.4) 

470  FORMATC  TEMP  Y  =  ',F10.4) 

480  FORMATC  TEMP  Z  =  ’,F10.4) 

450  FORMAT(IOX,I2,’  DXDYDZ',3F8.4) 

GO  TO  998 
C  POSTMORTEM 
999  CONTINUE 

WRITE(6,2999)NCOUNT,I,J,IX,IY,IZ,IP,ICI 
2999  FORMAT(/,"POSTMORTEM”,8I6,/) 

WRITE(6,*)RX(I),RY(I),RZ(I),VX(I),VY(I),VZ(I) 

GO  TO  3999 
998  CONTINUE 
STOP 
END 
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APPENDIX  II 


Reagents  procured  for  first  single  batch  (0.51@50%)  monodisperse  polystyrene  latex: 
Almog  dispersion  polymerization  method  (8),  prescription  according  to  Liu  et  al  (9,10) 


CHEMICAL  DESCRIPTION 

QUANTITY 

COST 

Poly(Ethylene  glycol)methyl  ether  (2000) 

250  g 

35.00 

1 1  -bromoundecanol 

50  g 

50.00 

p-toluenesulfonic  acid  monohydrate 

100  g 

6.00 

THE  anhydrous 

1  L 

37.00 

Ether 

2.5  L 

19.00 

Chloroform-distilled 

2.0  L 

17.00 

Methacryloyl  chloride 

100  ml 

24.00 

Azobisisobutyronitrile  (AIBN) 

100  g 

12.00 

Styrene 

2.5  L 

24.00 

Methanol 

2.5  L 

12.00 

Dichloromethane 

2.5  L 

15.00 

3 ,4-dihdro-2H-pyran 

100  ml 

10.00 

Magnesium  Sulphate 

500  g 

7.00 

Ethanol 

2.0  L 

20.00 

Potassium  Hydroxide 

Molecular  Sieves,  3A,  5A,  4A 

Triethylamine 

100  ml 

Toluene  anhydrous 

2.0  L 

35.00 

Chloroform  (CDCI3) 

2.0  L 

17.00 

